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Abstract 

This paper introduces the Inverse Gamma (IGa) stochastic volatility model with time-dependent 
parameters, defined by the volatility dynamics dVt = Kt {9t — Vt) dt + XtVtdBt- 

This non-affine model is much more realistic than classical affine models like the Heston stochastic 
volatility model, even though both are as parsimonious (only four stochastic parameters). Indeed, 
it provides more realistic volatility distribution and volatility paths, which translate in practice into 
more robust calibration and better hedging accuracy, explaining its popularity among practitioners. 

In order to price vanilla options with IGa volatility, we propose a closed-form volatility-of-volatility 
expansion. Specifically, the price of a European put option with IGa volatility is approximated by a 
Black-Scholes price plus a weighted combination of Black-Scholes greeks, where the weights depend 
only on the four time-dependent parameters of the model. 

This closed-form pricing method allows for very fast pricing and calibration to market data. The 
overall quality of the approximation is very good, as shown by several calibration tests on real-world 
market data where expansion prices are compared favorably with Monte Garlo simulation results. 
This paper shows that the IGa model is as simple, more realistic, easier to implement and faster to 
calibrate than classical transform-based affine models. We therefore hope that the present work will 
foster further research on non-affine models like the Inverse Gamma stochastic volatility model, all 
the more so as this robust model is of great interest to the industry. 

Key words: stochastic volatility, Inverse Gamma, volatility expansion, closed-form pricing, log¬ 
normal, mean-reverting SABR 
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1 Introduction 


The banking industry, especially in equity and foreign exchange, is currently experiencing a shift away 
from affine stochastic volatility models such as the Heston model, and towards non-affine stochastic 
volatility models such as the Inverse Gamma model. Non-affine stochastic volatility models have been 
shown to produce more realistic volatility paths and volatility distributions, to capture more accurately 
the dynamics of the market implied volatility surfaces, and to produce more reliable calibrations, thus 
reducing the realized volatility of delta-hedging P&Ls. 

Up until now, the popularity of affine models in spite of their empirical inadequacy has been due to one 
thing: tractability. Indeed, affine models provide quasi closed-form formulas for vanilla option prices by 
transform methods, in contrast to non-affine models. The purpose of the present work is to resolve this 
issue by presenting a fast pricing method for a non-affine stochastic volatility model. More precisely, 
we develop a closed-form expansion for the price of vanilla options under the non-affine Inverse Gamma 
stochastic volatility model, defined by the volatility dynamics dVt = Kt {Ot — Vt) dt + XtVtdBt. 
Implementing this new closed-form expansion is straightforward, and pricing speed is instantaneous. In 
fact, the closed-form expansion approach is much easier and much faster than the transform methods used 
for affine models. Moreover, our method is designed to deal naturally with time-dependent parameters. 
This freedom for the term structure of the model parameters makes the calibration process much easier 
for various maturities. 

We illustrate the accuracy of this closed-form expansion method on several foreign exchange market 
data sets. The speed and accuracy of the method make it ideal for industry use. The parameters 
generated through the fast calibration procedure can be used to directly price and hedge options under 
IGa stochastic volatility, but can also be used in the calibration of more general local-stochastic volatility 
models. 

The paper is organised as follows: 

• Section 2 defines the Inverse Gamma stochastic volatility model, discusses similar models in the 
literature, and discusses the advantages the Inverse Gamma model has over other classical one- 
factor stochastic volatility models. 

• Section 3 provides a closed-form volatility-of-volatility expansion for the price of a European put 
option under Inverse Gamma stochastic volatility. Furthermore, we provide an algorithm to easily 
compute the expansion coefficients for piecewise constant parameters, leading to fast calibrations. 

• Section 4 provides several numerical tests of the method on foreign exchange market data (AU- 
D/USD, USD/JPY, USD/SGD). In each example, the Inverse Gamma model is calibrated to the 
whole implied volatility surface and the expansion prices are compared to Monte Garlo prices. 
This allows us to assess both calibration error and expansion error. 

• Section 5 summarizes the presented methodology and provides some of our plans for future work 
in this area. 

2 The Inverse Gamma Stochastic volatility model 

This section defines the Inverse Gamma stochastic volatility model and discusses its properties. In this 
paper we use notations specific to foreign exchange (namely domestic and foreign interest rates), but 
note that the model itself is not limited to foreign exchange applications and can of course be readily 
used for other markets (equity, fixed-income, etc.) 

2.1 Definition 

Denote St and Vt as an exchange rate and its instantaneous volatility at time t, and T the time horizon 
considered. The dynamics of the Inverse Gamma (IGa) stochastic volatility model with time-dependent 
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parameters is given by 


dSt = {rd{t) - rf{t))Stdt + VtStdWt 
dVt = Kt {St ~ bt) dt + XiVtdBt 


( 2 . 1 ) 


d{W,B)^ = ptdt 


where {Wt^ ^t)o<t<T ^ two-dimensional correlated Brownian motion, = (rd(t))Q<(< 7 ' is the domestic 
interest rate, and r/ = (?'/(^))o<t<T foreign interest rate. There are four deterministic parameters: 

• K = {iit)o<t<T mean reversion of the volatility to the level 0. 

• 9 = (0t)o<t<T mean reversion level of the volatility. 

• A = (At)Q<-«-y is the volatility of volatility. 

• p = {pt)o<:t<T is the correlation between the respective Brownian motions of the underlying S and 
its volatility V. 

When the parameters are kept constant, the volatility (2.1) is driven by an Inverse Gamma process, 
yielding an inverse gamma distribution for the stationary distribution of volatility (cf. Appendix A. 2). 
Thus, we denote this model as the Inverse Gamma stoehastic volatility model (IGa model in short) with 
time-dependent parameters^. 

2.2 The IGa model in the literature 

A few classes of stochastic volatility models proposed in the literature contain the IGa model with 
constant parameters as a particular case. To make comparisons simpler, we use the same notations 
for the parameters of each class { k , 9 , A, p) and we remove the drift term from the dynamics of the 
underlying. 

• The Power Arch (or PARGH) stochastic volatility model (Fornari and Mele [2001]), 


dSt = VtStdWt 
dVf = k{9- VI) dt + XV[dBt 


corresponds to the IGa stochastic volatility model when p = 1. Remark that p = 2 corresponds to 
the GARGH diffusion model (cf. Table 2.1). 

• The Double Log-Normal stochastic volatility model (Gatheral [2007, 2008] , Henry-Labordere [2009] ) 
with its two cointegrated variance factors: 


dSt = ^/VtStdWt 

dVt = K (Vt - Vt'j dt + XVtdB, 


( 2 . 2 ) 



with correlations between the Brownian motions W, B and B . Indeed, the variance formulation 
of the IGa model can be reformulated as follows: 


dSt = VvStdWt 

dVt = (2K9Vt - [2k - A^] Vt) dt + 2XVtdBt 



which is a particular case of Double Log-Normal stochastic volatility with 100% correlation between 
B and B'. 


^See also Appendix B. 
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• The A — SABR modeP (Henry-Labordere [2008] Chapter 6 ), also known as mean-reverting SABR: 

dSt = VtS^dWt 

dVt = k{ 6- Vt) dt + XVtdBt 


One can see that the IGa model corresponds to the case /3 = 1, which, along with the case 
P = 1/2, is often considered in practice (see for example Shiraya and Takahashi [2011] or Shiraya 
and Takahashi [2014]). 

• The Generalized Inverse Gamma (GIGa) stochastic volatility model of Ma and Serota [2014]: 


dSt = VtStdWt 

dVt = K - Rt) dt + XVtdBt 


The special case 7=1 corresponds to an IGa diffusion for the volatility. 

• Finally, the closest model to (2.1) in the literature is the so-called “Log-normal Beta stochastic 
volatility model” of Sepp [2014, 2015]: 


dSt = VtStdWt 

dVt = K{e-Vt)dt + pVtdWt + eVtdBt (2.3) 

d{W,B)^ = 0 


The model (2.3) is in fact equivalent to the Inverse Gamma model (2.1), as if W and B are 
correlated Brownian motions (with correlation p), then B = pW + \/l — p'^W^ where is 
another Brownian motion, independent from W. Therefore (2.3) is equivalent to (2.1) with P = Xp 
and e = X\Jl — p^. For example, the typical equity case /3 « — 1 and e « 1 mentioned in Sepp [2014, 
2015] corresponds to a volatility of volatility A = ^/2 « 1.41 and a correlation p = —1/-\/2 « —0.71. 


Other classes of stochastic volatility that contain the IGa model include dVt = k {9 — Vt) dt + XVtdBt 
(Jerbi [2011], IGa when 77 = 1) and the general dVt = [q (t) Vp — s (t) Vp] dt+l (t) VP^^dBt (Itkin [2013], 
IGa when q = k9, s = k, I = X, a = 0, b = 1, though they focus on the closed forms that can be derived 
when 0 = 1 and 6 = 27 -|- 1, which excludes the IGa model). 

This list of models does suggest that the IGa model (2.1) is a sensible and reliable basis to model 
volatility, but also that the most efficient way to parsimonously enrich the model is not clear yet. This 
question is left for future research, and the rest of the paper will focus on the IGa model (2.1). 


2.3 Other models 

Over time, many stochastic volatility models have been proposed in the literature. Table 2.1 recalls some 
classical one factor stochastic volatility models with mean-reversion and correlation between volatility 
and underlying^ with constant parameters. To make comparisons to the IGa model easier, both volatility 
and variance formulations are given. 


^which, with our notations, is more accurately described as a k — SABR model. 
{W, B)f. = pdt in all models from Table 2.1 
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Name 


Variance formulation 


Volatility formulation 

dSt = {rd - rf)Stdt + VtStdWt 


dSt = {rd - rf)Stdt + y/VtStdWt 


Schobel-Zhu^ 

Heston^ 

3/2-model^ 
Log-Normal^’® 

GARCH® 

Inverse 

Gamma 


II 

{9- 

Vt) 

dt XdB^ 


dVt= ( 

HiO 

2 

8 

1 

^Vt)dt + 

\dBt 

dVt= ( 

fVt 

- 

2 8 

Vt^)dt + 

^Vt^dBt 

dVt= ( 

k9+^ 

Vt - nVt log(Vt)^ 

dt+XVtdB 

dVt= ( 

k9 1 

2 Vf 

- 

K 

2 8 

vyt+i 

^-VtdBt 

II 

{9- 

Vt) 

dt + XVtdBt 



dVt= {X^+2Ke^/Vi-2KVt)dt + 2X^/VidBt 

dVt = K {9 — Vt) dt + X\/VtdBt 

dVt = K {eVt - V^) dt + XV^ dBt 

dVt= {2[Ke+X^]Vt - KVt\og{Vt))dt+2XVtdBt 

dVt = k{9 — Vt) dt + XVtdBt 

dVt= {2K6y/Vt - [2K-X^]Vt)dt+2XVtdBt 


Table 2.1: One factor stochastic volatility models 

Among the models listed in Table 2.1, the models of Schobel-Zhu and Heston are affine, which means 
that the Fourier transform of the log-price can be computed explicitly. Because of their tractability, 
affine models have received a lot of attention in the literature, at the expense of the non-affine stochastic 
volatility models. Unfortunately, empirical analyzes suggest that the dynamics of market volatilities is 
much better described by non-affine models. Let us illustrate this point by comparing the Heston model 
and the Inverse Gamma model. 

Figure 2.1 displays the stationary distribution of the volatility under the Heston and Inverse Gamma 
models, with same mean (0.30) and same standard deviation (0.08 on Figures 2.1a and 2.1b, 0.16 on 
Figures 2.1c and 2.Id, 0.24 on Figures 2.1e and 2. If), using the results from Appendix A. 

On the one hand, the volatility distribution in the affine models from Table 2.1 (Schobel-Zhu, Heston) 
has a short right tail^, while it has a more realistic long right tail in the non-affine models (3/2 model, 
Log-Normal, GARGH, and Inverse Gamma). Figures 2.1b, 2.Id and 2.If (in log-scale), illustrate this 
difference between Heston and Inverse Gamma. The right tail of the Heston volatility decreases more 
quickly than that of the IGa volatility, therefore there is always a volatility level upon which the Heston 
volatility falls forever below the IGa one. 

The left tail is also better described with non-affine models. For example, with the Heston model, the 
volatility can reach zero if the parameters (k, 9, A, p) do not satisfy the Feller condition (2k, 9 j)? > 1). 
The effect of this condition on the left tail can be seen from Figures 2.1a and 2.1b (2 k 0/A^ = 3.63) to 
Figures 2.1c and 2.Id (2 k 0/A^ = 0.96) to Figures 2.1e and 2.If (2«;0/A^ = 0.49). One can clearly see 
how the distribution piles up close to zero, to the point where zero becomes the most likely value for the 
volatility (Figure 2.1e). Unfortunately, the Feller condition is almost always violated in practice (Glark 
[2011], da Fonseca and Grasselli [2011], Ribeiro and Poulsen [2013], ...), which means that Figure 2.1e 
represents the normal behavior of the Heston model on real data. 


^Schobel and Zhu [1999] 

^Heston [1993] 

^ Lewis [2000] 

^Wiggins [1987] 

®A more natural definition is dSt = Vd ~ rj)Stdt -|- e^*StdWt with dVt = k( 9 — Vt) dt + XdBt- See also Appendix B. 
® Lewis [2000] 

^using Rojo [1996] ’s tail classification for example 
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(b) mean = 0.30, standard deviation = 0.08, log scale 




Figure 2.1: Volatility density 
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2.3.1 Empirical evidence 

Figure 2.1 strongly suggests that the IGa model 2.1, while being as parsimonious as Heston (four 
parameters k, 0, X, p), provides a much better description of realized and implied market volatility. 
Indeed, a number of empirical studies suggest that a non-affine stochastic volatility model of the type 
(2.1) compares favorably to other possible models, especially to affine models such as the Heston model. 
We can summarize these studies as listed below. 

In equity markets: 

• In Bouchaud and Potters [2003] (Chapter 7), an empirical analysis of the volatility of the S&P 500 
in the period 1990-2001 is performed. The distribution of the volatility is shown to be accurately 
fitted by two distributions: a log-normal distribution, and an inverse gamma distribution. Overall, 
the best fit is provided by the inverse gamma distribution, especially on the right tail of the 
volatility distribution (Figures 7.7 p.118 and 7.8 p.ll9). When discussing the shortcomings of the 
Heston stochastic volatility models, the authors explicitly say that “the empirical distribution of 
the volatility is closer to an inverse gamma distribution than to a gamma distribution” (p.l43). 

• Gander and Stephens [2007] test several possible distributions for the volatility (Tempered stable. 
Generalized Inverse Gaussian (including Gamma), Positive Hyperbolic, Inverse Gaussian and In¬ 
verse Gamma) on 14 stocks traded on the New York Stock Exchange. It is shown that the Inverse 
Gamma distribution provides the best fit for option pricing. 

• Gatheral [2008] shows that the Double Lognormal stochastic volatility model (which contains 
(2.1) as a special case) fits SPX and VIX options much better than Double Heston, with stable 
parameters. 

• Christoffersen et al. [2010] show, on S&P 500 returns, VIX options and OTM S&P 500 index 
option data, between 1996 and 2004, that the GARCH diffusion model (Inverse Gamma variance) 
significantly outperforms several other models including the Heston and the 3/2 models. (They only 

consider the volatility parameterization (not variance) dVt = dt+ 

which does not contain the Inverse Gamma model.) 

• Furthermore, Kaeck and Alexander [2012] show that allowing for non-affine dynamics (like (2.1)) 
is more important than the inclusion of jumps. In particular, augmenting an affine model such as 
the Heston model with jumps leads to a stochastic volatility model that is still significantly inferior 
to a more parsimonious non-affine model like GARGH without jumps. 

• Finally, Ma and Serota [2014] analyze the volatility of the S&P 100, S&P 500 and DJIA, as 
well as the VIX index between 1990 and 2014. They show that the Generalized Inverse Gamma 
distribution (which contain the Inverse Gamma distribution) fits volatility best. 

In fixed-income markets: 

• Fornari and Mele [2001] calibrate the Power Arch model to futures contracts on the Italian 10-year 
government bond, between 1991 to 1997. The power p is estimated on three different subsamples. 
The results, 0.86, 0.99 and 1.19, are very close to an IGa volatility {p = 1). In particular, the 
GARGH diffusion model {p = 2) is rejected. 

• Fornari and Mele [2006] then study a slightly different stochastic volatility model, of the type 
dvr = k{ 9 — Vf)dt -I- XV^^dBt, with two parameters p and rj. They fit these two parameters to 
weekly 3-month US Treasury bills rates between 1973 and 1995. Their estimates, p = 1.0326 and 
fj = 1.0014, are statistically indistinguishable from 1 (Inverse Gamma volatility). Remark that the 
volatility term for the asset S (which is an interest rate in there case) is of the form Vty/StdWt, 
which is different from (2.1). However, they state that a volatility of the form Vt [S'jI dWt, d > 0.5 
(which contains (2.1) for d = 1), would not dramatically change their empirical results. 
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Beyond better calibration and more realistic volatility distribution and volatility paths, a model such 
as 2.1 generates a more realistic dynamics for the implied volatility surface (Tataru and Fisher [2012]), 
which reduces the tracking volatility in hedging P&Ls (Sepp [2015]), and drastically improves portfolio 
allocation (Hansis [2010]). 


3 Closed-form expansion for fast option pricing 

In the previous section, a number of reasons have been documented in favor of a non-affine Inverse 
Gamma stochastic volatility model over affine models. However, the main reason why affine models are 
used in practice is not their realism but their tractability. Indeed their Fourier (or Laplace) transform 
is available in closed form, which makes pricing possible by inverse transform. Non-affine models such 
as the IGa model do not have a closed form solution for their Fourier transform, which makes them a 
priori less tractable. 

In Sepp [2014], an affine moment-matching approximation of the moment generation function for (2.1) 
is proposed, making pricing possible by inverse transform. The advantage of this approach is that jumps 
can be factored in during the matching as well. Here we propose a more straightforward approach, 
namely a closed-form volatility-of-volatility expansion for vanilla options prices. Gompared to Sepp 
[2014], the main advantages of our approach are that: 

• It is simpler. If an approximation is to be made, it is more straightforward and intuitive to 
approximate the price directly rather than a transform of the price. This also makes pricing and 
subsequent calibration much faster. 

• No moment matching is required. Firstly, matching moments with an affine model can create 
unexpected problems (recall Figure 2.1e). Secondly, for most stochastic volatility models, moments 
higher than 1 cease to exist for large maturities (Andersen and Piterbarg [2007]). In Sepp [2014], 
the author advises against going beyond a second-order approximation, as only 5 moments were 
shown to exist. Our approach does not have such limitation on the order of approximation. 

• Importantly, our approximation approach is naturally suitable for time-dependent parameters. 

In the following subsection 3.1, we provide a closed-form expansion for the price of a European put option 
under Inverse Gamma volatility (2.1) with time-dependent parameters (equation (3.1)). The method is 
based on the methodology developed in Benhamou et al. [2010], adapted to the Inverse Gamma model 
and extended to time-dependent k. In Benhamou et al. [2010], the closed-form volatility of volatility 
expansion methodology was applied to the Heston stochastic volatility with time-dependent parameters. 
It was shown to be very accurate, and much faster than Fourier methods. 

The coefficients of the closed-form expansion are explicitly given by time integrals of the (time-dependent) 
parameters (equation (3.5)). Importantly, any shape for the time evolution of the parameters can be 
handled. In practice though, piecewise-constant parameters can be seen as a good compromise between 
richness and tractability. Thus, we present generic recusions formulas for the coefficients (3.5) when the 
model parameters are piecewise constant. These generic recursions are easy to implement, and can be 
used for expansion coefficients to any order. 

Finally, in order to compare our prices to another method, and as transform methods are not available, 
we explain how to implement an efficient Monte Garlo scheme to price options under IGa volatility. 


3.1 Closed-form expansion 


Theorem 3.1. The second-order expansion for the price Pjca = PiGa {Sq, K,T,r(i,rf-, k, 9, X, p) of a 
European put option with Inverse Gamma volatility is explicitly given by 


2 QiJrl 1 Q2iJr2 

PiGa = PbS {xo, 'tpT) + ^ &2i,T o (2^0, V't) + £ (3.1) 

i=o y i=o y 


with 


Xo 

Vo,t 

'IpT 
0^0,T 
0-1,T 
02,T 

bo,T 
62,T 


log(^o 

rt 




e ■’0 ^ I Vo 

r-T 


'0 


vltdt 




e^^^^'^^PsXsvlsds J e-^o"^"^vo,idt 

efo f^zdz / 2 ptXtVo^tdt / e~ ^0 >^^^^vo^udu 

J S J t 

’ Z*^ z*^ 

e^o^^<i^p^XsvlJs I e^o^^'^^ptXtvltdt / 


4 / Vo / e ^0 




oi.r 


(3.2) 

(3.3) 

(3.4) 


(3.5) 


where PBs{x,y) = Pbs ix,y; K,T,rd,rf) is the Black-Scholes put price with spot and integrated 
variance y, 


rT ( 1 f Ke~ ^0 1 ^ 


-e^e~ P 


( 1 , f 1 


(3.6) 


and £ is the error term of the second-order expansion. 


Proof. The proof or expansion (3.1) is available in Appendix C. It is based on the proof of the Heston 
expansion in Benhamou et al. [2010], that we extended to non-constant k, and adapted to the IGa 
model. □ 


Remark 3.1. One can also easily obtain closed-form expansion for the Greeks by adapting the proof of 
Theorem 3.1 to it. More generally the methodology adopted in Theorem 3.1 can be adapted to any 
options that have closed-form Black-Scholes prices with time dependent parameters, for example barrier 
options (Lo et al. [2003], Rapisarda [2003]). 
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3.2 Recursions for expansion coefficients 

The above coefficients (3.5) are expressed for general deterministic parameters {Kt,0t, Xt, Pt)o<:t<T- 
this subsection, we propose an explicit recursive algorithm to compute the coefficients when the param¬ 
eters are piecewise constant. 

Let To = 0 < Ti < T 2 < • • • < Tn = T he a partition of [0, T], We now suppose that the parameters are 
constant on each interval: 


Xt , Pt} 


{Ki,9i, Xi, Pi) Vt G [Tj, Ti_|_i [ 


(3.7) 


3.2.1 Integral operator 

Firstly, we define recursively the following integral operator ^ 


nT 

= y efo" Vte[0,T] 


^(k^,Z„).....(«i.Zi) - / VtG[0,r] 

Using this notation, the coefficients from the price expansion 3.1 can be expressed as follows: 


'07’ 
OO.T 
0.1,T 
02,T 
^0,T 


,(0’<.) 


‘"0,T 


= W, 


(2K,A^t)o_ ),(-2 k,1) 


0,T 

[k,p\vI),{-k,vo,,) 

rp 

(K.pAt)^_ ),(0.2pAi>o,'«o,.) „ (k,pAi>o_ ),(k,pAd^_ ),(-2k,1) 

rp ^C<Jq p 


= Auj, 


{2k,\^Vq)),{—k,vo,,),{—k,vo,. 


0,T 


(3.8) 


(3.9) 


(3.10) 


3.2.2 Recursions 

Let /i,Z 2 ,^ 3 , ■ ■ • be deterministic functions. When these functions are piecewise constant, 

lk{t) '■= lk,i , t G [Ti,Ti+i[ , fc = 1,2,... , 

then the following integral operators at time Tj_|_i can be expressed as a function of the integral operators 
at time Tp 




(niK,ZipJl) _ (niK,hv^l) (ni,0,pi) 


OjTi+i 


_ V ’0,./ I 7 

— ^0,Ti 


+1 


( 3 , 11 ) 


iV, 


n2K,l2VQ^ ) .(niKjZiPg^ 




= u;, 


n2K,l2Vg'^) ,(niK,liVg^) 


0,Ti 


+ 


(n2K,l2vl^) ,^(ni,0,pi) 

^0,Ti 


I „n2+nij j .^('‘2,11 


(n2,0,p2),(ni,0,pi} 


+1 


( 3 . 12 ) 


^which corresponds to Def. 5.1 in Benhamou et al. [2010] extended to non-constant 
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+ 


(n3K,Z3'»o^),(«2K,Z2t'o^),(niK,h'«o^) 

^o,r,+i 


+ 

+ 


(7131^,hvQ^ ), (n 2 K,Z 2 ^^o^ ), {n\ k,IiVq^ ) 
‘^0,Ti 


{n 3 K,l 3 V^^),[n 2 K,l 2 vll) (ni,0,pi) 


UsK.ZsJ^oi^) „2+ni7 7 ,^(»^2,0,p2),(ni,0,pi) 

^ 0 ,Ti ^ 0 ,Ti ''2,i>-l,iyTi,Ti+3 


„n3+n2+ni7 

^o,ri '■ 


7 7 ,^(?l 3 . 0 ,P 3 ).(”- 2 , 0 ,P 2 ),(ni, 0 ,pi) 

3,iZ2,iZl,i9^Ti,ri+i 


( 3 . 13 ) 


rt 1 

and so on, where ni,n 2 ,ri 3 ,... and Pi,P 2 ,P 3 , ■ ■ ■ are integers, eo,t = eJo and for all Ti < t < T^+i, 
p is defined as follows: 


(ni,mi,pi) 


,Ank,mk,Pk),- ,ini,mi,pi) 

^t,T,+3 


e-l It^ ^pij^ ^ f 




+1 


rik /y. Kzdz ( \mk ,,Pfc (nj,_i,mj,_i,ps,_i),--- ,(ni,mi,pi) , 
e * 7 ^ 0 ,sV^t,Ti+i 


^UkUiATais)^ (s)”*'' Vq'^ ^(^-l,mk-i,Pk-i),-- ,(ni,mi,pi)^^ 


where 7 (s) = with ATj = T^+i — Ti, and ni,mi,pi,..., nk,mk,Pk are integers. 

It should be noted that when the parameters are piecewise constant (equation (3.7)), p can be computed 
explicitly by recursion. Define Aki = k^+i — ki, A9i = 0^+1 — 9i, AXi = A^+i — Xi and Api = — pi, 

and let t € [Ti,Ti+i\. Then 


VQ,t = 0i + (fo,Ti - 


9i) 


1 


and, using the definition of p, basic integration and integration by parts, the following recursions hold 


n ,-(?il,mi,pi-l) I / _ Q \ /o'.”!' 

^i^t,Ti+3 + 1^0,T/ t/i) T’t,Ti+i 


(ni-l,mi,pi-l) 


Pi>0 


(ni,mi,pi) 




niKi 


ni = 0 , Pi = 0 


Til 7 ^ 0 , mi = 0 , Pi = 0 


^niKiATi gTij Kj 

niKi 


(ni,mi-l,0) 

niK^ATi^t,Ti+i 


rii 7 ^ 0 , mi > 0 , Pi 


0 


and for every integer fc > 1 : 
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Pfc > 0 






+ - 




AT’, 


,^(fife-i,"Xfc-l,Pfc-i)r-- I ATi _(nfe-i,mfc+m)c-i+l,Pfe-i),--- ^ _n ^ _n 

+ mk+l^t,Ti+i rifc - U, Pfc - U 


i + l 




TL}^ Ki 




+1 


I 1 _(ni;+TH._i,mfc_i,pfc_i),-' 




Aj=ol[^) jl [ukKiAT, 

i_ sr'^k rnj^ ( -1 \^k-3 (uk 

iKi ^j=0 j\ \nkKiATi / rt,Ti 


71]^ Ki 


rtik-J 


M + l 

_1 \'^k-3 J^nk+nk-i,mt,_x+j,Pk-i), 

i+l 


i'fik-l,mk-i,Pk-i),- 


Hk 7 ^ 0 , mfc = 0 , Pfc = 0 


Uk ^ 0, ruk > 0, pk = 0 


All these equation are sufficient to compute for any integers (ni, mi,pi), • • • , {nk,mk,Pk), 

making it possible to implement the integral recursions (3.11), (3.12) and (3.13), which, in turn, make 
it possible to compute explicitly the expansion coefficients (3.10) when the parameters of the stochastic 
volatility model are piecewise constant. 


3.3 Monte Carlo 

An obvious alternative method for computing the price of a European put option under the IGa stochastic 
volatility model is the Monte Carlo method. Using the expression 


PiGa {xo,vo) = E 


Pbs (xo+f ptVtdBt - ^ 


{ptVtfdt, r {I-pl)V,^dt 

Jo I 


(cf. equation (C.4)), only the simulation of the volatility is needed. Then, one can take advantage of 
the strong solution of the Inverse Gamma diffusion (Zhao [2009]) 


Vt = ^Vg + J KgOgZsds 
where Z is a geometric Brownian motion 

dZt = i^Kt + Ztdt — XtZtdBt 

ie. ^ ^ 

Zt = exp f ’ 

to derive the following unconditionally stable discretization scheme (called “Pathwise Adapted Lineariza¬ 
tion” in Kahl and Jackel [2006]) 

^ri t— -|- 2 '^?,! 

1 - 

^i„+i ^ -T- 

where 0 = to < in < iN = T is a time discretization of the interval [0,T], with At„ := 

tn+i — in and Ai?„ := — Bt^. This scheme ensures in particular that the paths of V remain 

positive. 


^ At„ — Xt^ABn 
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4 Numerical experiments 

This section provides numerical tests of the fast pricing method proposed in this paper. The Inverse 
Gamma stochastic volatility model will be calibrated to market data using the pricing formula (3.1). 
Three foreign exchange data sets will be used, detailed in subsection 4.1. For each test case, we will 
provide the calibration errors and expansion errors for implied volatility surfaces (Subsection 4.2). 

4.1 Datasets 

Three full sets of foreign exchange market data are provided here for easy benchmarking. We provide 
the strikes used for the implied volatility surface, and the equivalent constant rates for each maturity 
(the constant rate r^giT) equivalent to a time-dependent rate r{t), 0 < t < T is defined by req(T) := 
T fo The corresponding market implied volatility surfaces will be provided in subsection 4.2 

(Tables 4.4, 4.5 and 4.6), along with the corresponding calibration errors and expansion errors. All the 
numbers quoted here have been rounded. 

4.1.1 Set 1: AUD/USD 17 June 2014 {So = 0.9335) 


Strikes Rates 


Mat 

10 Put 

25 Put 

ATM 

25 Gall 

10 Gall 

Mat 

foreign 

domestic 

IM 

0.9103 

0.9233 

0.9356 

0.9469 

0.9572 

IM 

2.80% 

0 .21% 

3M 

0.8906 

0.9168 

0.9401 

0.9605 

0.9795 

3M 

2.89% 

0.31% 

6 M 

0.8664 

0.9100 

0.9469 

0.9780 

1.0078 

6 M 

3.03% 

0.45% 

lY 

0.8322 

0.9027 

0.9609 

1.0096 

1.0580 

lY 

3.26% 

0.69% 


Table 4.1: Market data, AUDUSD, 17 June 2014 


4.1.2 Set 2: USD/JPY 11 June 2014 (S'o = 102.00) 


Strikes Rates 


Mat 

10 Put 

25 Put 

ATM 

25 Gall 

10 Gall 

Mat 

foreign 

domestic 

IM 

99.78 

100.88 

101.99 

103.09 

104.16 

IM 

0 .20% 

-0.04% 

3M 

97.47 

99.77 

101.98 

104.15 

106.31 

3M 

0.29% 

0.07% 

6 M 

94.75 

98.47 

102.00 

105.46 

109.06 

6 M 

0.40% 

0.16% 

lY 

90.04 

96.34 

102.01 

107.67 

114.06 

lY 

0.52% 

0 .21% 


Table 4.2: Market data, USDJPY, 11 June 2014 

4.1.3 Set 3: USD/SGD 04 September 2014 (^o = 1.2541) 


Strikes Rates 


Mat 

10 Put 

25 Put 

ATM 

25 Gall 

10 Gall 

Mat 

foreign 

domestic 

IM 

1.2397 

1.2466 

1.2542 

1.2637 

1.2755 

IM 

0.16% 

0.17% 

2M 

1.2334 

1.2432 

1.2542 

1.2688 

1.2871 

2M 

0.19% 

0.19% 

3M 

1.2286 

1.2406 

1.2543 

1.2729 

1.2970 

3M 

0.28% 

0.27% 

6 M 

1.2152 

1.2339 

1.2545 

1.2836 

1.3233 

6 M 

0.48% 

0.47% 

lY 

1.1945 

1.2232 

1.2548 

1.3018 

1.3704 

lY 

0.61% 

0.57% 


Table 4.3: Market data, USDSGD, 04 September 2014 
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Our empirical analyses suggest that these three data sets are representative of the behavior on market 
data of the model and of the expansion. In particular, our empirical findings will be similar for each 
data set, leading us to present only three examples, as providing more examples would not provide much 
additional information. 


4.2 Calibration 

Using the closed-form expansion (3.1), we can respectively calibrate the Inverse Gamma model to the 
three market data sets. The four (piecewise-constant) stochastic parameters k, 9, A and p are calibrated, 
along with the initial volatility Vq- This calibration process can be assessed from the implied volatility 
calibration error. Then, using these parameters, we will estimate the expansion error, by comparing the 
implied volatility provided by the expansion (3.1) to a Monte Carlo price (Subsection 3.3), computed 
with 24 time steps per day and M = 1 000 000 paths (to keep both bias and variance very low). 

4.2.1 Set 1: AUD/USD 17 June 2014 {Sq = 0.9335) 

The estimated initial volatility is Vq = 6.49%, the estimated piecewise-constant stochastic parameters 
are given by 



K 

9 

A 

P 

IM 

4.19 

6.39% 

1.71 

-0.40 

3M 

2.33 

11 . 01 % 

1.12 

-0.74 , 

6 M 

2.26 

11.85% 

1.25 

-0.73 

lY 

1.80 

12.52% 

0.87 

-0.92 


and the calibration and expansion errors are given below (rounded to the nearest basis point): 



10 Put 

25 Put 

ATM 


25 Call 

10 Call 

IM 

7.48 [- 0 . 04 ] [- 0 . 04 ] 

6.87 [-0.02] [-0.01] 

6.38 [ 

0 . 06 ] [ 0 . 00 ] 

6.19 [ 0.00] [ 0.02] 

6.19 

[- 0 . 03 ] [- 0 . 00 ] 

3M 

8.46 [ 0 . 05 ] [- 0 . 12 ] 

7.48 [- 0 . 05 ] [- 0 . 03 ] 

6.68 [ 

0 . 03 ] [ 0 . 03 ] 

6.36 [- 0 . 04 ] [ 0 . 08 ] 

6.39 

[ 0 . 05 ] [- 0 . 08 ] 

6 M 

9.93 [ 0 . 02 ] [- 0 . 15 ] 

8.43 [- 0 . 03 ] [- 0 . 03 ] 

7.30 [ 

0 . 06 ] [ 0 . 04 ] 

6.82 [- 0 . 07 ] [ 0 . 13 ] 

6.90 

[ 0 . 08 ] [- 0 . 29 ] 

lY 

11.51 [- 0 . 05 ] [- 0 . 19 ] 

9.53 [ 0 . 00 ] [- 0 . 04 ] 

8.05 [ 

0 . 16 ] [ 0 . 07 ] 

7.47 [- 0 . 14 ] [ 0 . 19 ] 

7.57 

[ 0 . 15 ] [- 0 . 63 ] 


Table 4.4: AUDUSD Market implied volatility [calibration error] [expansion error] in % 


The median absolute deviation of the calibration error is 5.0bp, and its mean absolute deviation is 5.7bp. 
The median absolute deviation of the expansion error is 5.5bp, and its mean absolute deviation is 10.9bp. 
Overall, the median absolute deviation of the total error is 6.0bp, and its mean absolute deviation is 
10.5bp. 

4.2.2 Set 2: USD/JPY 11 June 2014 {Sq = 102.00) 

The estimated initial volatility is Vq = 4.42%, the estimated piecewise-constant stochastic parameters 
are given by 
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K, 

e 

A 

P 

IM 

8.23 

7.96% 

2.47 

- 0.10 

3M 

5.00 

6.47% 

1.32 

-0.19 

6 M 

3.62 

9.32% 

1.61 

-0.15 

lY 

2.10 

6.74% 

1.88 

- 0.22 


and the calibration and expansion errors are given below (rounded to the nearest basis point): 



10 Put 

25 Put 

ATM 


25 Call 


10 Call 

IM 

6.08 

[ 0.08] [-0.00] 

5.76 

[ 0.00] [-0.01] 

5.53 [ 

0.01] [ 0.00] 

5.63 

[-0.07] [ 

0.00] 

5.81 

[-0.02] [-0.00] 

3M 

7.08 

[ 0.02] [-0.02] 

6.53 

[-0.03] [-0.03] 

6.15 [ 

0.01] [ 0.00] 

6.20 

[-0.03] [ 

0.01] 

6.44 

[ 0.04] [-0.01] 

6 M 

8.19 

[-0.12] [ 0.06] 

7.42 

[ 0.02] [-0.03] 

6.95 [ 

0.10] [-0.01] 

7.00 

[ 0.02] [ 

0.00] 

7.36 

[-0.05] [ 0.02] 

lY 

9.90 

[-0.15] [ 0.32] 

8.61 

[ 0.13] [-0.12] 

7.95 [ 

0.07] [-0.07] 

8.04 

[-0.04] [- 

0.02] 

8.69 

[-0.06] [ 0.08] 


Table 4.5: USDJPY Market implied volatility [calibration error] [expansion error] in % 

The median absolute deviation of the calibration error is 4.0bp, and its mean absolute deviation is 5.4bp. 
The median absolute deviation of the expansion error is 1.5bp, and its mean absolute deviation is 4.1bp. 
Overall, the median absolute deviation of the total error is 2.0bp, and its mean absolute deviation is 
3.9bp. 

4.2.3 Set 3: USD/SGD 04 September 2014 (^o = 1.2541) 

The estimated initial volatility is Vb = 3.16%, the estimated piecewise-constant stochastic parameters 
are given by 



tv 

6 

A 

P 

IM 

2.90 

4.03% 

2.30 

0.49 

2M 

2.88 

4.19% 

1.64 

0.49 

3M 

2.85 

4.44% 

2.37 

0.58 

6 M 

2.76 

4.08% 

1.68 

0.51 

lY 

2.81 

4.27% 

2.31 

0.67 


and the calibration and expansion errors are given below (rounded to the nearest basis point): 


10 Put 25 Put ATM 25 Call 10 Call 


IM 

3.06 

[-0.07] [ 

0.05] 

3.02 1 

: 0.06] [ 

0.02] 

3.30 

[ 0.04] [ 0.01] 

3.79 

[ 0.01] [-0.03] 

4.43 

[-0.05] [-0.02] 

2M 

3.15 

[-0.02] [ 

0.07] 

3.15 1 

: 0.02] [ 

0.05] 

3.45 

[ 0.01] [-0.00] 

4.10 

[-0.02] [-0.07] 

4.84 

[ 0.00] [-0.03] 

3M 

3.23 

[-0.00] [ 

0.07] 

3.24 1 

: 0.00] [ 

0.06] 

3.57 

[ 0.02] [-0.00] 

4.35 

[-0.02] [-0.09] 

5.20 

[-0.01] [ 0.02] 

6 M 

3.52 

[-0.03] [ 

0.08] 

3.47 1 

-0.02] [ 

0.09] 

3.80 

[ 0.05] [ 0.02] 

4.78 

[ 0.02] [-0.11] 

5.85 

[-0.07] [ 0.17] 

lY 

3.84 

[-0.05] [ 

0.04] 

3.78 1 

-0.04] [ 

0.11] 

4.20 

[ 0.12] [ 0.03] 

5.40 

[ 0.15] [-0.07] 

6.78 

[-0.19] [ 0.68] 


Table 4.6: USDSGD Market implied volatility [calibration error] [expansion error] in % 


The median absolute deviation of the calibration error is 2.0bp, and its mean absolute deviation is 4.4bp. 
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The median absolute deviation of the expansion error is 5.0bp, and its mean absolute deviation is S.Obp. 
Overall, the median absolute deviation of the total error is T.Obp, and its mean absolute deviation is 
7.8bp. 

4.3 Comments 

As shown by Tables 4.4, 4.5 and 4.6, both calibration errors and expansion errors are small overall (only 
a few basis points of absolute deviation). 

The overall calibrations are very good. The worst cases occur for lY maturity and 25 Call strike on Set 
1 (+19bp), for lY maturity and 10 Put strike on Set 2 (—15bp), and for lY maturity and 10 Call strike 
on Set 3 (—19bp). In other words, they occur for long maturities and strikes far out of the money. 

As to the expansion error, it is as expected very small for short maturities and close to the money 
(ATM), but can get larger for long maturities and far from the money. The worst cases occur for lY 
maturity and 10 Call strike on Set 1 (—63bp), for lY maturity and 10 Put strike on Set 2 (+32bp), and 
again for lY maturity and 10 Call strike on Set 3 (+68bp). Though these worst case errors can look 
large in implied volatility terms, the absolute option price values are very small far out of the money, 
meaning that the actual error is not that large in absolute price terms. Overall, we observe that the 
expansion error increases with maturity T and volatility of volatility A (as expected from Benhamou 
et al. [2010]), and with absolute correlation jpj (compare Set 2 to the two other sets). 

To decrease the calibration errors, one could try more general stochastic volatility models, as discussed in 
Subsection 2.2, though adding too many parameters may generate overfitting and damage the stability 
and robustness of the model. 

To decrease the expansion error, the most obvious solution is to compute a higher-order expansion 
(expansion (3.1) is only a second-order expansion). Another idea is to take advantage of the fact that 
the expansion error is much smaller in practice when the stochastic parameters are constant (only a 
few bps everywhere). Benhamou et al. [2010] tried to exploit this idea by first calibrating a model with 
constant parameters for each maturity, and then turning these calibrations into an “equivalent” model 
with piecewise constant parameters. However, the equivalent piecewise constant model is determined 
using the same second-order expansion, therefore this two-step calibration procedure does not improve 
the expansion error of the final model with piece wise-constant parameters. 

Depending on the intended application, the accuracy of the present second-order expansion (3.1), as 
observed in Subsection 4.2, may be sufficient. For example, if one only needs the calibrated stochastic 
parameters to feed into a more general local-stochastic Inverse Gamma volatility model, then the local 
volatility component can easily eliminate the residual price discrepancies (cf. Sepp [2014]). 

4.4 Comparison to Heston 

Finally, we compare the Inverse Gamma calibration to a classical Heston stochastic volatility calibration. 
Subsection 2.3 provided theoretical reasons to favor the Inverse Gamma model over the Heston model, so 
it is interesting to check if, for example, one of the two models provides significantly better calibrations 
in practice, as this is one of the important criteria for practical use by the industry. We use the same 
three datasets, and calibrate the Heston model using the semi-closed-form of Heston [1993]. 

4.4.1 Set 1: AUD/USD 17 June 2014 {Sq = 0.9335) 

For the Heston model, the estimated initial variance is Vq = 0.41%, the estimated piecewise constant 
stochastic parameters are given by 
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n 

9 

A 

P 

2kQI\^ 

IM 

1.16 

1.28% 

0.32 

-0.32 

0.30 

3M 

0.97 

2.32% 

0.48 

-0.49 

0.20 , 

6 M 

1.01 

1 .88% 

0.51 

-0.54 

0.15 

lY 

1.02 

1.85% 

0.51 

-0.52 

0.14 


and the calibration error is given below (rounded to the nearest basis point): 


10 Put 25 Put ATM 25 Call 10 Call 


IM 

7.48 

[ 0.02] 

6.87 

[-0.03] 

6.38 

[ 0.00] 

6.19 

[-0.04] 

6.19 

[-0.05] 

3M 

8.46 

[ 0.03] 

7.48 

[-0.07] 

6.68 

[ 0.01] 

6.36 

[ 0.03] 

6.39 

[ 0.07] 

6 M 

9.93 

[ 0.03] 

8.43 

[-0.05] 

7.30 

[ 0.03] 

6.82 

[ 0.02] 

6.90 

[-0.06] 

lY 

11.51 

[-0.04] 

9.53 

[-0.10] 

8.05 

[ 0.10] 

7.47 

[ 0.08] 

7.57 

[-0.07] 


Table 4.7: AUDUSD Market implied volatility [Heston calibration error] in %, 


The median absolute deviation of the calibration error is 4.3bp, and its mean absolute deviation is 4.6bp. 


4.4.2 Set 2: USD/JPY 11 June 2014 (^o = 102.00) 

The estimated initial variance is Vg = 0.28%, the estimated piecewise constant stochastic parameters 
are given by 



K 

9 

A 

P 

2k9I\^ 

IM 

1.17 

1.39% 

0.23 

-0.11 

0.62 

3M 

1.10 

1.76% 

0.40 

-0.21 

0.25 . 

6 M 

1.09 

1.73% 

0.44 

-0.22 

0.20 

lY 

1.04 

1.92% 

0.48 

-0.43 

0.17 


and the calibration error is given below (rounded to the nearest basis point): 



10 Put 

25 Put 

ATM 

25 Call 

10 Call 

IM 

6.08 [-0.06] 

5.76 [-0.06] 

5.53 [ 0.00] 

5.63 [-0.10] 

5.81 

[-0.13] 

3M 

7.08 [-0.00] 

6.53 [-0.06] 

6.15 [ 0.00] 

6.20 [-0.02] 

6.44 

[ 0.04] 

6 M 

8.19 [ 0.13] 

7.42 [-0.04] 

6.95 [-0.05] 

7.00 [-0.05] 

7.36 

[ 0.06] 

lY 

9.90 [ 0.04] 

8.61 [-0.03] 

7.95 [-0.02] 

8.04 [-0.01] 

8.69 

[ 0.02] 


Table 4.8: USDJPY Market implied volatility [Heston calibration error] in % 


The median absolute deviation of the calibration error is 4.0bp, and its mean absolute deviation is 4.6bp. 
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4.4.3 Set 3: USD/SGD 04 September 2014 (^o = 1.2541) 


The estimated initial variance is Vq = 0.11%, the estimated piecewise constant stochastic parameters 
are given by 



K 

0 

A 

P 

2 k0 /A2 

IM 

1.25 

0.49% 

0.28 

0.41 

0.16 

2M 

1.15 

0.62% 

0.49 

0.29 

0.06 

3M 

1.24 

0.54% 

0.48 

0.23 

0.06 

6 M 

1.51 

0.32% 

0.52 

0.44 

0.04 

lY 

1.21 

0 .86% 

0.49 

-0.26 

0.09 


and the calibration error is given below (rounded to the nearest basis point): 


10 Put 25 Put ATM 25 Call 10 Call 


IM 

3.06 

[-0.07] 

3.02 

[ 0.02] 

3.30 

[ 0.00] 

3.79 

[-0.02] 

4.43 

[-0.02] 

2M 

3.15 

[-0.05] 

3.15 

[-0.01] 

3.45 

[-0.01] 

4.10 

[-0.05] 

4.84 

[ 0.01] 

3M 

3.23 

[ 0.02] 

3.24 

[-0.00] 

3.57 

[-0.01] 

4.35 

[-0.10] 

5.20 

[-0.00] 

6 M 

3.52 

[ 0.01] 

3.47 

[-0.03] 

3.80 

[ 0.03] 

4.78 

[-0.04] 

5.85 

[ 0.13] 

lY 

3.84 

[ 0.09] 

3.78 

[-0.04] 

4.20 

[-0.02] 

5.40 

[-0.08] 

6.78 

[ 0.13] 


Table 4.9: USDSGD Market implied volatility [Heston calibration error] in % 


The median absolute deviation of the calibration error is 2.5bp, and its mean absolute deviation is 4.0bp. 
4.4.4 Result analysis and comparison 

For the three examples, the calibrated parameters do not satisfy the Feller condition: the Feller ratio is 
much smaller than 1. This is the reason why we could not rely on the Heston closed-form expansion of 
Benhamou et al. [2010] for calibration. Indeed, when the Feller ratio is so low, the expansion error can 
be massive. For the USD/JPY parameters, the average absolute expansion error is around 80bp. For the 
AUD/USD and USD/SGD parameters, the error is so large that some expansion prices can even become 
negative. When the implied volatility exists, the average absolute expansion error is around 200bp 
for both AUD/USD and USD/SGD parameters. These large errors were expected, as the expansion 
approach for the Heston model is not expected to work when the Feller condition is not satisfied (the 
error analysis in Benhamou et al. [2010] requires the Feller condition to hold). This means that the 
expansion scheme for the Heston model cannot be used in practice because, as recalled in subsection 
2.3, the Feller condition is virtually never satisfied on real-world market data. For this reason, we used 
the slower semi-closed-form pricing formula of Heston for calibrating the model. 

The table below summarizes the calibration error (median and mean, in basis points) for both Inverse 
Gamma and Heston stochastic volatility models. 


AUD/USD USD/JPY USD/SGD 


calibration error (bp) 

median 

mean 

median 

mean 

median 

mean 

Inverse Gamma 

5.0 

5.7 

4.0 

5.4 

2.0 

4.4 

Heston 

4.3 

4.6 

4.0 

4.6 

2.5 

4.0 
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One can see that the calibration error is of the same order. If anything, the Heston calibration is slightly 
better (—0.8bp on average). Therefore, the calibration error itself does not provide clear and sufficient 
indication on which stochastic volatility model is a better model in practice. To assess the quality of 
a stochastic volatility model, the stability of the estimated parameters (the “variance” in statistical 
terms) is as relevant as the calibration error (the “bias’” in statistical terms). Indeed, a very flexible 
stochastic volatility model with many parameters is likely to have a very small calibration error on a 
given market implied volatility curve. However, the parameters of such a model are also likely to be very 
unstable over time, if recalibrated every day as is common in practice. The stability of parameters is an 
important criterion for dynamic hedging to work in practice and for the option price to be meaningful at 
all. This “bias-variance tradeoff” discussion suggests that more numerical tests, in particular stability 
tests, are needed to further compare the Inverse Gamma and the Heston stochastic volatility models. In 
particular, a dynamic hedging backtest using the two models for several options over a long time period 
would be useful. We plan to carry out such tests in the future. 

At this stage of the comparison, we recommend the Inverse Gamma model over the Heston model for 
foreign exchange option pricing for at least three reasons: 

• Firstly, even though, based on the three examples studied in this paper, the Heston calibration 
seems slightly more accurate, the fact that the optimal parameters are very far from satisfying the 
Feller condition is a major issue, as the very unrealistic volatility distribution resulting from this 
breach (see Figure 2.1e) is likely to affect the quality of dynamic hedging and the pricing of more 
exotic options. 

• Secondly, from our experience of using the model for more than forty currency pairs, the closed- 
form expansion for Inverse Gamma vanilla options (Theorem 3.1) works quite well on real-world 
market data (unlike the Heston expansion). If needed, the accuracy of the expansion, already 
good, can always be enhanced using a higher-order expansion (e.g. 4th order). 

• Lastly, the closed-form expansion approach is the fastest method for calibration purposes, and 
therefore is highly desirable in practice when calibration needs to be performed for hundreds of 
currency pairs or securities on at least a daily basis. 

5 Conclusion 

This paper has introduced the Inverse Gamma stochastic volatility model, as defined by the volatility 
dynamics dVt = Kt {Ot — Vt) dt + XtVtdBt. 

The volatility distribution in this model is more consistent with market dynamics than alternative one 
factor affine stochastic volatility models such as the Heston model. 

We have proposed a closed-form volatility of volatility expansion for the price of a European put option 
under this stochastic volatility model, and simple and straightforward recursion formulae to instanta¬ 
neously compute the coefficients of the expansion when the four time-dependent stochastic parameters 
K, 6, A and p are piecewise constant. 

We have demonstrated the viability of the second-order expansion scheme on three test cases from foreign 
exchange (AUD/USD, USD/JPY and USD/SGD). Both calibration error and expansion error are small 
overall (only a few basis points of absolute deviation on average). 

Potential improvements can be made and have also been discussed in the paper, such as introducing 
additional stochastic parameters, computing higher order terms in the expansion, performing more com¬ 
prehensive calibration backtests, and studying the local-stochastic volatility version of this model. We 
hope to foster further academic research on these non-affine models favoured by industry practition¬ 
ers, as this paper demonstrates that they not only provide a more accurate representation of market 
dynamics, but also still be tractable thanks to expansion methods. 
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A Stationary distribution of volatility 


A.l Heston 


For the Heston stochastic volatility model with constant coefficients, 

dSt = {rd-rf)Stdt+y^StdWt (A.l) 

dVt = K {9 — Vt) dt + X\/VtdBt 

d{W,B)^ = pdt, (A.2) 


the long-run distribution of the variance Vt is a Gamma distribution with parameters kr = fdO and 
0r = ^, where /3 := Remark that the probability density function of Vt can reach 0 if /30 < 1 (Feller 
condition). 

Consequently, the long-run Heston volatility Yt = \/Vt has a generalized Chi distribution 


Px v) 


1 





with parameters a = 0, b = and v = 2(59. 

In particular, the moments of the long-run volatility Yt = x/Vt are given by: 


E[Yt] 

E 


r(fcr + ^) 

F(fcr) 


\/^ 


krOr = 0 


F m x/f39 


(A.3) 


A.2 Inverse Gamma 

For the Inverse Gamma stochastic volatility model with constant coefficients, 

dSt = {vd - rf)Stdt + VtStdWt 
dVt = k{9- Vt) dt + XVtdBt 
d{W,B)^ = pdt, 

the long-run distribution of the volatility process Vt is an inverse Gamma distribution 

I /3“ -a-l 

Pr-i = :p7-wa; e - 

F(a) 

with parameters aj—i = 1 + 13 and /^r-i = (d9, where j3 := "^ (see Barone-Adesi et al. [2005]). Therefore 


E[Rt] 

Var [y] 


Pt- 


CUp-i — 1 


= 9 {iff 13 >0) 




(op-i — 1) (op-i — 2) (3 1 


(*/// 3 >l) 


(A.4) 


B The Log-Normal terminology 

There exists three different volatility dynamics that have been called log-normal in the literature: 


dVt 

= nVtdt + XVtdBt 

(B.l) 

dVt 

= K {9 -Vt)dt + XVtdBt 

(B.2) 

dlogiVt) 

= K {9 - log{Vt)) dt + XdBt 

(B.3) 
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In the first formulation (B.l), the volatility is modeled by a geometric Brownian motion, which is 
log-normally distributed. It is a special case of SABR model (with (3=1). It is not mean-reverting. 

In the last formulation (B.3), the logarithm of the volatility is an Ornstein-Uhlenbeck process, which 
has a normally-distributed stationary distribution. Therefore the long-term distribution of the volatility 
Vt is indeed log-normally distributed. 

The intermediate formulation (B.2) combines characteristics from a geometric Brownian motion (the 
volatility of volatility Wt is proportional to the volatility Vt) and from an exponential Ornstein- 
Uhlenbeck (the mean-reversion effect towards a level 6). However the stationary distribution of the 
volatility is not a log-normal distribution but an inverse Gamma distribution (Barone-Adesi et al. [2005], 
Zhao [2009], Sepp [2014, 2015]). 

Practitioners sometimes refer to (B.2) as a “Log-Normal” or “mean-reverting Log-Normal” stochastic 
volatility model, but as this terminology can be ambiguous and misleading, we choose to call it Inverse 
Gamma stochastic volatility model, which is consistent with the stationary distribution of the volatility 
(B.2). 


C Proof of main expansion 


In the IGa model (2.1), factoring out the drift rates rj, and ry, the dynamics of the log-spot Xt reads 

dXf = — ^dt+VtdWt, Xq = xq 

dVt = Kt{9t-Vt)dt + \tVtdBt, Uo = t'o (C.l) 

d (VP, B)^ = ptdt , 


where xq = log (^o). Define a perturbed process (A®, V^) as follows 

dXl = -^^dt + Vt^dWt, X^ = xo 

dU/ = Kt{0t-Vt^)dt + eXtV^^dBt, Vo^ = vo 
d{W,B)^ = Ptdt, 


(G.2) 


and define 


g{e) = e:xjp J rd(<)dt^ E - exp J {rd{t) - r f{t)) d^ 


■Al 


4-1 


(G.3) 


so that g(l) = PjGat the price of European put with IGa stochastic volatility, which is the quantity we 
want to compute. Remark that g(0) reduces to a Black-Scholes price. 

The expression (G.3) for g{e) can be simplified. Remark that Wt can be decomposed into 


Wt = ptBt + ^l-pjdBi 
where B^ is a Brownian motion independent from B. Therefore 

X! = xo + PtVt^dBt - i f\yn^dt + 1^ \/^- piVt^dBi . 

Let JF^ = Qct<T filtration generated by B. One can see that A|. has a Gaussian distri¬ 

bution with mean xq -V ptV^dBt - ^ J^(Vf)'^dt = |xo -V ptVfdBt - ^ (ptVf)'^dt'j - ^ /^^(l - 
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(it and variance (1 — Pf)(V/)^dt. Consequently, 


5 (e) = exp(-J rd{t)dtjE E (x-expl-j {rd{t) - rf{t)) dtj + X^j \F^ 


= E 


/ T T T N 

PbsUo + I^ PtV,^dBt-^l^ {ptVn^dt,!^ {l-p^,){vn^dt 


(C.4) 


where PBsi^i u) is the Black-Scholes price of a European put option with spot and integrated variance 
y (equation (3.6)). 


For any non-negative integer n, define := and Vnp ■= 


. One can check that 


e :=0 


and 


dVod = Kt{6t - VQ t)dt + eXtVo fdBt, Vq q = vq 

dV^,t = -i^tV^,tdi + nXtV^_iddBt + eXtV^fdBt, Vl^ = 0 , n > 1 


= e- do ^ Ksdsedo 

Vn,t = e-do-^d. f edS^^d.nX^Vn-i,sdB,, n>l 
Jo 


(C.5) 

(C. 6 ) 


where vq^ ■= Vq^ is deterministic. 

At this point, the main idea of the proof is to approximate 5 ( 1 ) using a Taylor expansion of g{e) around 
e = 0, as 5 and its derivatives reduce to a Black-Scholes formula when e = 0. 

Thus, we approximate Vt = V^^ (e = 1) using = vo,t (^ = 0) by applying the Taylor formula to the 


function £ 1 —>■ 


Vt = vo,t + ki,t -I- 2^2,4 + . •. 


Doing the same to the function e 1 —>■ (E/)^ yields 

(Vtf = {vo,tf + 2vo,tVpt + {vo,tV2,t + (Vptf) + ... 
To simplify notations, define, for * > 0 and j > 0, 


d^+^pBS d^+^PBS 

dx'^y^ dx^yi 


^0 + Jo PtVo,tdBt Pi'^ltdt ’ “ Pl)vl,tdtj ■ (C.7) 

Then, the second-order Taylor expansion of 5 (e) around e = 0, valued at e = 1, reads (keeping only 
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second-order terms) 


5(1) = E[Pbs] 


(C. 8 ) 


(a :=) 

{Cy ■■=) 

(a, :=) 

{Cyy :=) 

iC.y ■■=) 


E 


E 






-h E 
+ £ 


dPss 

dx 

dpBS 

dy 

d^pBS 

dx^ 

d^pBS 

dy^ 

d^Pps 

dxdy 


J Pt {^i,t + ~ 2 f + 'Uo,tV2^t) dt 

j (l ~ Pt) (2'Co,tl4,t + + 'l' 0 ,tV 2 ,i) dtj' 

J Pt {Vip dBt ~ 2 J Pt (^■I'o.tl^i.t) dt 


{1- p'^t){2vo,tViPdt\ 


J PtiViPdBt-^J p'} {2vo,tViP dt\ I J {1 - p'f) {2vo,tViP dt 


where £ is the expansion error. From definition (C.7), taking e = 0 in (C.4) shows that 


E [Pbs] = Bps ^^o,y fo^td^ 


Similarly, 


E 


Id^+^pBs] 

d^+^Pps / 

dx'^yl 

dx'^yl \ 


^0, J vltdt\ . 


(C.9) 


(C.IO) 


The next part of this appendix is devoted to the computation of all the constants C^, Cy, Cxx, Cyy and 
Cxy Our main tools will be the following relation between Black-Sholes greeks: 


dPps 

dy 


{x,y) = - 


1 fd^Pps 


dPps 




the following product identity between two stochastic processes X and Y : 

rT 


XtYt = XoFo + [ {XtdYt + YtdXt + d (X, Y),) , 
^0 


(C.ll) 


(C.12) 


and the following Lemma: 

Lemma C.l. If C = I PuVo^udBu^ for a differentiable function I with derivative , then its first 

Malliavin derivative (G) = {Df (G))^^^ is equal to Df (G) = l^^'^ ^ Jq PuVo,udBu^ ps'Co,sl{s < T}. 
Therefore, using Lemma 5.2 from Benhamou et al. [2010], 


E 



PutJo,udB. 



= E 




Puttt),udB. 



Pst^O.stXsds 


(C.13) 


for 0 < t < T. Most of the time, we will use (C.13) with t = T. 


• Computation of Cx- 
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Using (C.13), 

dpRS ^ 


E 


dx 


Pt Vi,t + -U2.t dB, 


= E 


d^pBS 

dx^ 


Therefore, using (C.ll), 

dpBS ^ 


= E 


dy 


Pt (2vo,tVi^t + vo,tV2P dt 


-E 


Pt {vo,tyi,t + dt 


dpBS \ ^ 2 xr 2 


dx I 2 


pM,tdt 


• Computation of Cxx- 
Using (C.12), then (C.13), and finally (C.ll), 

Cxx = E 


O^Pbs 

dx^ 


lo Uo 



= E 

= 2E 


d^pBS d^pBS^ 


dx^ dx"^ 
\ d^pBS [^ ( P 

dxdy Jo \Jo 
Therefore, so far, using (C.ll), 


PsVi^sdBs - ! plvo,sVi^sds^ [ptVx^tdBt - plvo,tVi^tdt] 
PsVi^sdBs - j plvo,sVi,sds^ p1vo,tVi,tdt 


-E 

2 


d^pBS 2..2 ,, 



0 lio 


-E 

2 


VPbS P 2,,2 ,, 


/*^ 1 1 3 ^P Q c'^ 

PsVi^sdBs-J plvo^sVi^sdsj PtVo^tVi^tdt +-E J PtVltdi 


Cx + Cy + Cxx - E 

{Cxy-=) + 2E 

• Computation of Cxy'- 
Using (C.12), 

d^P, 


dpBS 

dy 


d^pBS 

dxdy Jo Uo 


i^vojVi^t + Upt + t'o,tV2,t) dt 

PsVi^sdBs - J plvo,sVi^sds^ pfvojVi^tdt 


■T ( nt 



Cxy — E 




PtVi^tdBt-J^ p?z;o.tUi,tdt| (1 - p?) ( 2 z;o,tUi.t) dt 

PsVi^sdBs - J p'^gVo,sVi^sds^ (l - Pt) (2vo,tUi,t) dt 



~dxiy Jo {/ ~ P^s) {‘^vo,sVi^s) ds^ {ptVi^tdBt - PtVojVi^tdt} 


Using (C.13), 


-'xy 



PsVi^sdBs- J plvo,sVi^sds^ {2vojVi^t)dt 

(l - pf) {2vo,sVi^s) dsj {ptVi^tdBt - PtVojVi.tdt} 
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Define G = vo^tVi^t- Then Df (G), its first Malliavin derivative w.r.t. B, is given by 

D^iG) = + (,.,,) 

^ fn Hzdz K^dz \ -fi f ^ 

^ ^ vo,te -^0 ^ e-^o ^ AsZ;o,sl{s < 


= vodVid ^^ PsVo^s'^is <T}-\- vo^te 

ox^oy oxoy 


using the definition of Vi^t (C.6). Therefore, using Lemma 5.2 in Benhamou et al. [2010], 
^d^pBS 

-fo.tn.t ( ^ 

dt 


2 / E 
10 
pT 

= 2 

'o 


pT r f^ 2 p / rt 

L ^ hsf "-‘'M (/„ 


dt 


= 2 


^ E g|^ PsVi^sdB, 

[ E^/ psVi,sDf (G) ds 
Jo Jo 


dt 


= 2E 


= 2E 


o^Pbs r f r 2 


dx'^dy JO v^o 
'^BS 


dx'^dy 


Ps'^o,sVi,sdsj VojVijdt 

j plvo^sVi^sds^ vojVijdt 


+ 2 r Vo, te-p[ 
Jo I dxoy Jq 


fp. p pi 

/ eP-^‘i-X,p,vo,sV,,,ds 
dxoy Jq 


dt 


gop rT / rt 


-E 


OPbs 

dy 


vo,tV2,tdt 


where we used, for the last equality, equation (C.13) and the definition of V 2 ,t (equation (C.6)). There¬ 
fore, using (C.13) and then (C.ll), 


Gxy + Cxy — 4E 


4E 


-h E 


d^pBS 

r( 

1- 

to 

Jo V 

'd^pBs 

n 

_1 


Ps'>Jo,sVi,sds vo,tVi,tdt 


dpBS 

dy 


vo,tV2,tdt 


Using (C.12), 


Thus 


Gyy = 4E 


O^Pbs ^ 

dy"^ 


J |y - pI) vo,sVi,sds^ {1 - pf) Vo, tVi,tdt 


Gxy Cxy Gyy - 4E 


d^Pi 


BS 


= 2E 


d^pBS 

dy^ 


( / vo^sVi.sds I vodVid^t 

-hE 

\Jo J 



vo,tVi,tdt 


-E 


QPbs 

dy 


dPss f 
dy Jo 


vo,tV2,tdt 


vo,tV2,tdt 


At this stage, combining all the terms together, the second-order expansion for g{l) (equation (C.8)) 
simply becomes 


g{l) = Pbs (xo, J Vo,tdt\ +E {2vo,tVi,t + + 2 vo,tV 2 ,t) dt 


2E 


d^pBS 

dy'^ 


vo,tVi,tdt ^ 
(C.14) 
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Equation (C.14) is much simpler than equation (C.8), but still contains the stochastic processes 
and E 2 ,t- The last part of the proof is to simplify (C.14) further to get an expansion with explicit, 
deterministic coefficients. 

Recall the short-hand notation for deterministic integrals from equations (3.8) and (3.9): 


UJ. 


'^t,T 


(K'n4n),....(Klll) 

t,T 


= e^o" vte[o,r] 

. 


Vt G [0, T] . 


We will use extensively the following Lemma: 

Lemma C.2. (Lemma 5.4 in Benhamou et al. [2010]) For any deterministic integrable function f and 
any continuous semimartingale Z such that Zq = 0, 

[ f{t)Ztdt= [ uff^UZt 
Jo Jo 

Proof. Apply Ito’s lemma to the product Zt. □ 

The computation of each type of expectation in equation (C.14) is summarized in the following Lemma. 
Lemma C.3. The following equalities hold 


E 


E 


E 


I ^J ptVo^tdBtj j PtVi^tdt 

I ( [ PtVo.tdBt] [ l3tV2,tdt 


I ( J PtVo.tdBt j J PtVltdt 


{k,p\v^),{-k,P) 
‘^O.T 


E 


l^^'‘ ( J Pt'^o,tdBt 


= w, 


{K,pXvg),(0,2p\vo,.),{ — K,fi) 
0,T 

{2k,^X^Vq ),( — 2 k .,/ 3 ) 


E 


( / PtVo,tdBt 


(C.15) 


(C.16) 


•^0,T 


E 


I / PtVo,tdBt 


2^('^>P^’'o,.).(K,pAt)o_4,(-2K./3)g 


l^^'> f J ptVo.tdBt 


(C.17) 


E 


I / PtVo.tdBt < / PtVi^tdt 




I / PtVo.tdBt 




{K,pXVg ),{-K,P) 
‘^O.T 


E 


( / PtVo.tdBt 


(C.18) 


where fd, a deterministic function, and I, a twice-differentiable function, are such that these expectations 
exist. 

Proof. Using Lemma C.2 and equation (C.13), 


E 


I ( [ PtVo,tdBt] [ PtVi,tdt 


= E 


= E 


l(^J^ PtVo,tdBt'^ edo-^'i^X„vo,sdBsdt 

l(^l^ PtVo,tdBt'^ 

1 [ PtVo,tdBt 


{k,pXvI^_),{-k,/ 3 ) 

^0,T 


28 






























E 


^ ( [ PtVo,tdBt\ [ l3tV2,tdt 


= E 


= E 


PtVo,tdB}j Pte-^0 

PtVo,tdBt] f J-T’^^e^o-^dz 2 p^X^vo,tVxtdt 


= E 


= E 


'E 


l\^j PtVo,tdBtj j J-T’^^e^o-‘dz2x^Vi,tdB, 

ptVo^tdB^ j uj[~r^'^^ 2 pt\tVQ^t J e^o ’^^'^^XsVo^sdBsdt 

(^l\tVo,tdB}j , ,(«.pA<,),(o,2pA.„,.).(-../3) 

Using (C.12), 

PtVo.tdBt'^ dt 

\ e^o '^^<i^X,vo,,dB,dt 


( J Pt'^o,tdBt 


E 


I ( J ptVo^tdBt j J PtVltdt 


= E 


= E 


+ E 


I / PtVo,tdBt / Pte 


lij ptVo^tdBtj j Pte ^-^0 



First, using Lemma C.2, 

rT \ rT 


E 


\/o I e2/o«^<i-A2uo%ds|dt 


(2k,A^Do_ ),(-2k,/3) 

^0,T 


E 


I I J ptVo,tdBt 


Then, using Lemma C.2, Lemma (C.l) and then equation (C.15), 

rt 


E 


/N PtVo,tdBt\ j 2|y 


= E 


= E 


l(^J^ PtVo,tdBt'^ 244 ''’'^^e/o-^'^-Atz;o,t|^ e^o "*‘'"A,uo.sdB,| dS* 
( r PtVo,tdBt] r do-^d^p^XtvltV.^tdt 


= 2u}, 


T \ .T 

PtVo^tdBt j 

0 / ^0 

{ k , p\vq ),(/t,pAi;Q ),( —2«;,/3) 


0,T 


E 




yielding (C.17). Finally, using (C.12) and Lemma C.2, 


E 



^ ( f PtVo,tdBt] < f PtVi^tdt 
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Using equation (C.17), 


2E 


l( [ ptVo,tdBt] [ 




I / PtVo,tdBt 




( / PtVo,tdBt 


Then, using Lemma C.l, Lemma C.2, and equation (C.15), 
2E 


l\ I ptVo,tdBt] f l3sVi,sds'^edo-^^^XtVo,tdBt 


= 2E 


= 2E 


PsVi^sds^dt 

( r PtVo,tdB}\ 




( J Pt'^o,tdBt 


(K,pAVg ),(K,pAVg_),(-K,P),(-K,P) „ («,P''''«0,. )■(-'=._ J (k ,p\Vg ) ,(-K ,P) 


+ 2W| 


= S W, 


^0,T 


and using that 4 wq y 
(Benhamou et al. [2010] p. 34) yields equation (C.18). 

To conclude the proof, apply Lemma C.3 and equation (C.IO) to each term of equation (C.14), yielding 


ff(l) — Pbs ^2;o, j +E J (2uo,tVi^t + + 2uo,tV2_t) (it 

Bbs ^^0; J' 


+ 2E 


Q^Pbs 

dy"^ 


vo.tVi,tdt 


w, 


(k.pAdq ),(-k,2«o,,) d^PsS 
dxdy 


0,T 


Xo, / Vo ,dt 


''0,T 


a^o, / Vo,t 


{2k,,Pvo ),{-2k, 1) dPsS 
dy 


(k,pAi>0_ ),(k,pAuo ),(-2k,1) 3^PbS 
dx'^dy 


dt 


+ ° 


Xo, / vltdt 




(k.pAuq ),(0,2pAi;o,.),(-K,2wo,.) d^PsS 

dx'^dy 


0,T 


Xo, / Vofdt 


(2K,A^j;g ),(-K,vo,.),i-K,vo,.) O'^PbS 

4‘^o,t 


+ ^ 1 


dy 

(k,pAko_ ),(-K,uo,.) 1 Q'^PbS 


^ a;o, / Vo^t 


^ dt 


dx'^dy'^ > jQ 

which corresponds to the announced Theorem 3.1. 


Xo, / vl.dt 
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